/******************************************************************
 *      Equation Parser Utility Class                             *
 *                                                                *
 *  Please report bugs, incorrect coding, enhancements and        *
 *  updates to the wizardscripts forum                            *
 *  at www.wizardscripts.com/forum                                *
 *                                                                *
 *	Rights: The author retains all rights to this script      *
 ******************************************************************/

/** 
 * @main    Plik z definicja klasy Matrix
 * @file    macierz.cpp
 * @author  www.wizardscripts.com/forum
 * @modified Grzegorz Sejnota
 *
 * @created on 19 marca 2009, 10:41
 * @version 0.0001
 */


#include "macierz.h"
#include "wyjatki.h"
#include <iostream>
#include <cmath>
#include <deque>
#include <numeric>
using namespace std;

//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//++++++++++++++ POCZATEK [definicji klasy Matrix] +++++++++++++++++++++++++++++
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

//======= konstruktor bezargumentowy ===========================================
/**
 * Konstruktor domyslny, bezargumentowy. Tworzy macierz o wymiarach 1x1 o wartosci
 * 0.0.
 */
Macierz::Macierz() {
    macierz = vector< vector<double> > (1, vector<double> (1, 0));
}

//======= konstruktor ==========================================================
/**
 * Konstruktor klasy z vectora vectorow. Tworzy macierz na podstawie podanego 
 * vectora vectorow, wypelniajac pola liczby kolumn i wierszy.
 *      @param m - vector vectorow wartosciami macierzy
 */
Macierz::Macierz(const vector< vector<double> > &m)
    : macierz(m) {
}

//======= konstruktor konwertujacy z typu double ===============================
/**
 * Konstruktor konwertujacy z typu double.
 *      @param m - vector vectorow wartosciami macierzy
 */
Macierz::Macierz(double liczba){
    macierz = vector< vector<double> > (1, vector<double> (1, liczba));
}

//======= konstruktor macierzy o podanych wymiarach ============================
/**
 * Konstruktor macierzy o podanych wymiarach i wypelniajacy macierz podana wartoscia.
 *      @param wie - liczba wierszy
 *      @param kol - liczba kolumn
 *      @param liczba - wypelnienie
 */
Macierz::Macierz(int wie, int kol, double liczba){
    macierz = vector< vector<double> > (wie, vector<double> (kol, liczba));
}

//======= konstruktor ==========================================================
/**
 * Konstruktor klasy z vector. Tworzy macierz (wektor) na podstawie podanego 
 * vectora. W zaleznosci od argumentu flagi wierszowy utworzony wektor jest wektorem
 * wierszowym lub kolumnowym (domyslnie wierszowy)
 *      @param w - vector z wartosciami wektora
 *      @param wierszowy - flaga oznaczajaca czy tworzony wektor ma byc wierszowy
 *                  czy kolumnowy?
 */
Macierz::Macierz(const vector<double> &w, bool wierszowy){
    if (wierszowy){
        macierz = vector< vector<double> > (1, w);
    }else{
        macierz = vector< vector<double> > (w.size(), vector<double>(1, 0));
        for(int i=0; i<liczbaWierszy(); i++){
            macierz.at(i).at(0) = w.at(i);
        }
    }
}

//======= konstruktor z deque ==================================================
/**
 * Konstruktor klasy z deque. Tworzy macierz (wektor) na podstawie podanej kolejki. 
 * W zaleznosci od argumentu flagi wierszowy utworzony wektor jest wektorem
 * wierszowym lub kolumnowym (domyslnie wierszowy)
 *      @param w - deque z wartosciami wektora
 *      @param wierszowy - flaga oznaczajaca czy tworzony wektor ma byc wierszowy
 *                  czy kolumnowy?
 */
Macierz::Macierz(const deque<double> &d, bool wierszowy){
    if (wierszowy){
        macierz = vector< vector<double> > (1, vector<double>(d.size(), 0.0));
        for(int i=0; i<liczbaKolumn(); i++){
            macierz.at(0).at(i) = d.at(i);
        }
    }else{
        macierz = vector< vector<double> > (d.size(), vector<double>(1, 0.0));
        for(int i=0; i<liczbaWierszy(); i++){
            macierz.at(i).at(0) = d.at(i);
        }
    }
}

//======= konstruktor kopiujacy ================================================
/**
 * Konstruktor kopiujacy klasy Macierz
 *      @param &wzor - referencja do obiektu klasy Matrix, 
 *                                        ktory ma byc skopiowany.
 */
Macierz::Macierz(const Macierz & wzor) {
    this->macierz = wzor.macierz;
}

//======= destruktor  ==========================================================
/**
 * Destruktor klasy Matrix.
 */
Macierz::~Macierz() {
}

//======= metoda statyczna I ==================================================
/**
 * Metoda statyczna zwracajaca macierz jednostkowa o podanym rozmiarze.
 *      @param rozmiar - rozmiar macierzy jednostkowej
 *      @return macierz jednostkowa o rozmiarze rozmiar
 */
Macierz Macierz::I(int rozmiar){
    Macierz m(rozmiar, rozmiar, 0);
    for(int i = 0; i<rozmiar; i++){
        m.at(i, i) = 1.0;
    }
    return m;
}

//======= przeladowany operator przypisania ====================================
/** 
 * Przeladowany operator przypisania =. 
 *
 *     @param wzor - macierz ktora ma byc przypisana
 *     @return referencja do macierzy wynikowej (do ktorej przypisano wzor)
 */
Macierz & Macierz::operator=(const Macierz &wzor) {
    this->macierz = wzor.macierz;
    return *this;
}

//======= przeladowany operator wpisania do strumienia =========================
/** 
 * Funkcja przeładowanego operatora wypisanie wielomianu do strumienia
 *
 *     @param output_stream - referencja do strumienia wyjściowego 
 *     @param print_matrix - macierz do wypisania  
 *     @return output_stream - referencja do strumienia
 */
ostream & operator<<(ostream& output_stream, const Macierz& print_matrix) {
    for (int r = 0; r < print_matrix.liczbaWierszy(); ++r) {
        for (int c = 0; c < print_matrix.liczbaKolumn(); ++c) {
            output_stream << print_matrix(r, c) << " ";
        }
        output_stream << endl;
    }
    return output_stream;
}

//======= przeladowany operator () =============================================
/** 
 * Funkcja przeładowanego operatora funkcyjnego. Operator zwraca element macierz
 * na pozycji podanej w argumentach. W wypadku odwolania sie do nieistniejacego
 * elementu zwraca wyjatek. Zwracana wartosc stanowi L-wartosc. @see at().
 *
 *     @param row - numer wiersza (indeksowany od 0)
 *     @param col - numer kolumny (indeksowany od 0) 
 *     @return wartosc macierzy na pozycju row col
 */
double Macierz::operator() (int row, int col) const{
    if (row >= liczbaWierszy() || col >= liczbaKolumn() )
        throw Wyjatek("Proba odwolania sie do nieistniejacego elementu");
    return this->macierz[row][col];
}

//======= przeladowany operator [] =============================================
/** 
 * Funkcja przeładowanego operatora tablicowego. Operator zwraca referencje do wektora 
 * o indeksie podanym w argumencie. Operator umozliwa wywolania Macierz m[][].
 * Nie posiada jednak kontroli zakresu, w przeciwienstwie do metody at(int, int)
 * @see at(int, int).
 *     @param ind - numer wiersza (indeksowany od 0)
 *     @return referencja do wiersza(wektora) o indeksie ind
 */
vector<double> & Macierz::operator[] (int ind){
    return this->macierz[ind];
}

//======= przeładowany operator += =============================================
/**  
 *  Dodaje do macierzy macierz m i zwraca wskaźnik do macierzy wynikowej.
 *      @param m - macierz, który ma zostać dodany
 *      @return macierz wynikowowa (*this + m)
 */
Macierz Macierz::operator+=(const Macierz & m) {
    *this = *this +m;
    return *this;
}

//======= przeładowany operator -= =============================================
/**  
 *  Odejmuje do macierzy macierz m i zwraca wskaźnik do macierzy wynikowej
 *  wynikowego
 *      @param m - macierz, ktora ma zostac odjeta
 *      @return macierz wynikowowa (*this - m)
 */
Macierz Macierz::operator-=(const Macierz & m) {
    *this = *this -m;
    return *this;
}

//======= przeładowany operator *= =============================================
/**  
 *  Mnozy macierz this przez macierz m i zwraca wskaźnik do macierzy wynikowej
 *  wynikowego
 *      @param m - macierz z prawej strony mnozenia
 *      @return macierz wynikowowa (*this * m)
 */
Macierz Macierz::operator*=(const Macierz & m) {
    *this = *this * m;
    return *this;
}

//======= przeładowany operator /= =============================================
/**  
 *  Dzieli macierz this przez macierz m i zwraca wskaźnik do macierzy wynikowej
 *      @param m - macierz z prawej strony mnozenia
 *      @return macierz wynikowowa (*this * inv(m) )
 */
Macierz Macierz::operator/=(const Macierz & m) {
    *this = *this / m;
    return *this;
}

//======= funkcja liczbaWierszy ================================================
/**  
 * Zwraca liczbe wierszy macierzy.
 *      @return liczba wierszy macierzy
 */
int Macierz::liczbaWierszy() const {
    return this->macierz.size();
}

//======= funkcja liczbaKolumn =================================================
/**  
 * Zwraca liczbe kolumn macierzy.
 *      @return liczba kolumn macierzy
 */
int Macierz::liczbaKolumn() const {
    return this->macierz.at(0).size();
}

//======= funkcja obliczWyznacznik =============================================
/**  
 * Oblicza wyznacznik macierzy metoda eliminacji Gaussa.
 *      @return wyznacznik
 */
double Macierz::obliczWyznacznik() {
    Macierz matrix(this->macierz);
    Macierz new_matrix(this->macierz);
    Macierz temp_matrix(this->macierz);

    int SIZE = 0;

    if (matrix.liczbaWierszy() != matrix.liczbaKolumn()) {
        throw Wyjatek("Macierz nie jest kwadratowa, nie mozna obliczyc wyznacznika");
    } else {
        SIZE = matrix.liczbaWierszy();
    }
    
    if (SIZE == 1){
        //cout<<"MAMY SKLAR!"<<endl;
        //cout<<"wyznacznik = "<<new_matrix.at(0, 0);
        return new_matrix.at(0, 0);
    }

    if (SIZE == 2) {
        return ( (matrix(0, 0) * matrix(1, 1)) - (matrix(0, 1) * matrix(1, 0)));
    } else {
        int pivot_row = 0;
        int pivot_col = 0;
        double pivot_point = matrix(pivot_row, pivot_col);

        bool last_pivot = false;
        bool first_pivot = true;
        while (!last_pivot) {
            for (int row = (pivot_row + 1); row < SIZE; ++row) {
                for (int col = 0; col < SIZE; ++col) {
                    //Check if last pivot.
                    if (pivot_point == 0) {
                        new_matrix.at(row, col) = 0;
                    } else {
                        if (first_pivot != false) {

                            new_matrix.at(row, col) = matrix(row, col) - matrix(row, pivot_col) *
                                    matrix(pivot_row, col) / pivot_point;
                        } else {

                            new_matrix.at(row, col) = temp_matrix(row, col) - temp_matrix(row, pivot_col) *
                                    temp_matrix(pivot_row, col) / pivot_point;
                        }
                    }
                }
            }
            if (pivot_row == (SIZE - 2)) {
                last_pivot = true;
                break;
            }
            //Move the pivot up one
            pivot_row = pivot_row + 1;
            pivot_col = pivot_col + 1;
            pivot_point = new_matrix(pivot_row, pivot_col);
            first_pivot = false;
            temp_matrix = new_matrix;
        }

        //compute the determinant from the new matrix
        double determinant = 1;
        for (int r = 0; r < SIZE; ++r) {
            determinant = determinant * new_matrix(r, r);
        }

        return determinant;
    }
    return 0;
}

//======= funkcja obliczSlad ===================================================
/**  
 * Oblicza slad macierzy. Slad moze zostac obliczony dla macierzy kwadratowych, 
 * jezeli macierz nie jest kwadratowa funkcja wyrzuca wyjatek.
 *      @return slad
 */
double Macierz::obliczSlad() {
    if (liczbaKolumn() != liczbaWierszy())
        throw Wyjatek("Nie mozna obliczyc sladu, macierz nie jest kwadratowa");
    double slad = 0.0;
    for(int w = 0; w<liczbaWierszy(); w++){
        slad += macierz.at(w).at(w);
    }
    return slad;
}

//======= funkcja odwroc =======================================================
/** 
 * Funkcja odwracajaca macierz. Funkcja wykorzystuje metode adjoint.
 * W przypadku gdy macierz jest nieodwracalna lub 
 * nie mozna jej odwrocic bo nie jest kwadratowa funkcja zwraca wyjatek.
 *     @return macierz odwrotna
 */
Macierz Macierz::odwroc() const {
    int old_rows = this->liczbaWierszy();
    int old_cols = this->liczbaKolumn();

    vector< vector<double> > vect_matrix(old_rows, vector<double>(old_cols, 0));
    Macierz new_matrix(vect_matrix);
    Macierz temp(this->macierz);

    //if the determinant is zero or the matrix is not square,
    //-then it does not have an inverse. 
    if ((abs(temp.obliczWyznacznik()) <= 0.0001) || (old_rows != old_cols)) {
        throw Wyjatek("Nie mozna odwrocic macierzy");
    }
    
    if (old_rows == 1 && old_cols == 1){
        //cout<<"MAMY SKLAR!"<<endl;
        new_matrix.at(0, 0) = 1/temp.at(0, 0);
        return new_matrix;
    }

    //check if two x two square matrix
    if (new_matrix.liczbaWierszy() == 2) {
        double temp_zero = temp(0, 0);
        double temp_one = temp(0, 1);
        double determinant = temp.obliczWyznacznik();
        //change left diagnol
        new_matrix.at(0, 0) = temp(1, 1);
        new_matrix.at(1, 1) = temp_zero;

        //change right diagnol
        new_matrix.at(0, 1) = temp(1, 0);
        new_matrix.at(1, 0) = temp_one;

        for (int rows = 0; rows < 2; ++rows) {
            for (int cols = 0; cols < 2; ++cols) {
                if( (rows+cols)%2 == 0 ){
                    new_matrix.at(rows, cols) = new_matrix(rows, cols) / determinant;
                }else{
                    new_matrix.at(rows, cols) = -new_matrix(rows, cols) / determinant;
                }
            }
        }

    } else {
        //if matrix is larger than 2x2 use the adjoint method
        vector< vector<double> > cofactor_matrix(old_rows - 1, vector<double>(old_cols - 1, 0));
        Macierz cofactor(cofactor_matrix);
        double determinant = temp.obliczWyznacznik();
        //get cofactor matrix
        //-a cofactor matrix is the original matrix excluding
        //-the current elements row AND column.
        for (int cof_row = 0; cof_row < old_rows; ++cof_row) {
            for (int cof_col = 0; cof_col < old_cols; ++cof_col) {

                for (int r = 0; r < old_rows; ++r) {
                    for (int c = 0; c < old_cols; ++c) {
                        if ((cof_row != r) && (cof_col != c)) {
                            if ((r < cof_row) && (c < cof_col)) {
                                cofactor.at(r, c) = temp(r, c);
                            } else if ((r < cof_row) && (c > cof_col)) {
                                cofactor.at(r, c - 1) = temp(r, c);
                            } else if ((r > cof_row) && (c < cof_col)) {
                                cofactor.at(r - 1, c) = temp(r, c);
                            } else {
                                cofactor.at(r - 1, c - 1) = temp(r, c);
                            }
                        }
                    }
                }
                //the new value in the matrix to be returned is the
                //-determinant of the cofactor matrix.  +/- signs are
                //-used interchangably
                //cout<<"indeks"<<cof_row<<cof_col<<cofactor<<cofactor.obliczWyznacznik()<<endl;
                if( (cof_row+cof_col)%2 == 0 ){
                    new_matrix.at(cof_row, cof_col) = (cofactor.obliczWyznacznik())/determinant;
                    //cout<<"wynik"<<cof_row<<cof_col<<" = "<<new_matrix.at(cof_row, cof_col)<<endl;
                }else{
                    new_matrix.at(cof_row, cof_col) = -(cofactor.obliczWyznacznik())/determinant;
                    //cout<<"wynik"<<cof_row<<cof_col<<" = "<<new_matrix.at(cof_row, cof_col)<<endl;
                }
            }
        }
    }
    //before returning the matrix, it needs to be transposed
    return (new_matrix.transponuj());
}


//======= funkcja transponuj ===================================================
/**  
 * Transponuje macierz.
 *      @return macierz transponowana
 */
Macierz Macierz::transponuj() const {
    Macierz original_matrix(this->macierz);
    Macierz new_matrix(liczbaKolumn(), liczbaWierszy(), 0.0);

    for (int r = 0; r < new_matrix.liczbaWierszy(); ++r) {
        for (int c = 0; c < new_matrix.liczbaKolumn(); ++c) {
            new_matrix.at(r, c) = original_matrix(c, r);
        }
    }
    return new_matrix;
}

//======= funkcja at() =========================================================
/** 
 * Funkcja pobierajaca referencje do elementu macierzy na pozycji podanej w argumencie. 
 * W wypadku odwolania sie do nieistniejacego elementu zwraca wyjatek. 
 * Zwracana wartosc stanowi L-wartosci. @see operator().
 *     @param row - numer wiersza (indeksowany od 0)
 *     @param col - numer kolumny (indeksowany od 0) 
 *     @return wartosc macierzy na pozycju row col
 */
double &Macierz::at(int row, int col) {
    if (row >= liczbaWierszy() || col >= liczbaKolumn())
        throw Wyjatek("Proba odwolania sie do nieistniejacego elementu");
    return macierz[row][col];
}

//======= funkcja wiersz() =====================================================
/** 
 * Funkcja pobierajaca wiersz o indeksie podanym w argumencie. Pobrany wiersz
 * nie moze stanowic L-wartosci.
 *     @param ind - numer wiersza (indeksowany od 0)
 *     @return wiersz o indeksie ind
 */
vector<double> Macierz::wiersz(int ind) const{
    if (ind >= liczbaWierszy())
        throw Wyjatek("Proba odwolania sie do nieistniejacego elementu");
    return this->macierz[ind];
}

//======= funkcja kolumna() ====================================================
/** 
 * Funkcja pobierajaca kolumne o indeksie podanym w argumencie. Pobrana kolumna
 * nie moze stanowic L-wartosci.
 *     @param ind - numer kolumny (indeksowany od 0)
 *     @return wiersz o indeksie ind
 */
vector<double> Macierz::kolumna(int ind) const{
    if (ind >= liczbaKolumn())
        throw Wyjatek("Proba odwolania sie do nieistniejacego elementu");
    return this->transponuj().wiersz(ind);
}
//------------------------------------------------------------------------------
//-------------- KONIEC [definicji klasy Matrix] -------------------------------
//------------------------------------------------------------------------------



//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//++++++++++++++ POCZĄTEK [deklaracji funkcji pomocniczych] ++++++++++++++++++++
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

//======= operator == ===========================================================
/** 
 * Funkcja przeładowanego operatora porownania.
 *     @param m1 - macierz pierwsza
 *     @param m2 - macierz druga
 *     @return <code>true</code> jesli macierze sa sobie rowne
 *             <code>false</code> w przeciwnym razie
 */
bool operator==(const Macierz &m1, const Macierz &m2) {
    // ----- macierze maja rozny rozmiar ---------------------------------------
    if (m1.liczbaKolumn() != m2.liczbaKolumn() || m1.liczbaWierszy() != m2.liczbaWierszy()) {
        return false;
    }
    // ----- sprawdzenie rownosci elementow ------------------------------------
    for (int w = 0; w < m1.liczbaWierszy(); w++) {
        for (int k = 0; k < m1.liczbaKolumn(); k++) {
            if (m1(w, k) != m2(w, k)) return false; //nie sa rowne, false
        }
    }
    return true;
}

//======= operator != ==========================================================
/** 
 * Funkcja przeładowanego operatora porownania.
 *     @param m1 - macierz pierwsza
 *     @param m2 - macierz druga
 *     @return <code>true</code> jesli macierze nie sa sobie rowne
 *             <code>false</code> w przeciwnym razie
 */
bool operator!=(const Macierz &m1, const Macierz &m2) {
    // ----- macierze maja rozny rozmiar ---------------------------------------
    if (m1.liczbaKolumn() != m2.liczbaKolumn() || m1.liczbaWierszy() != m2.liczbaWierszy()) {
        return true;
    }
    // ----- sprawdzenie rownosci elementow ------------------------------------
    for (int w = 0; w < m1.liczbaWierszy(); w++) {
        for (int k = 0; k < m1.liczbaKolumn(); k++) {
            if (m1(w, k) != m2(w, k)) return true; //nie sa rowne, false
        }
    }
    return false;
}


//======= operator + ===========================================================
/** 
 * Funkcja przeładowanego operatora dodawania macierzy. Jezeli jedna z macierzy
 * jest skalarem to zostanie wykonane dodwanie skalara do kazdego elementu 
 * drugiej macierzy. Jezeli zadna z macierzy nie jest skalarem i nie zgadzaja
 * sie ich wymiary funkcja wyrzuca wyjatek klasy Wyjatek.
 *     @param w1  -  składnik pierwszy 
 *     @param w2  -  składnik drugi
 *     @return true - jesli wielomiany sa sobie rowne
 */
Macierz operator+ ( const Macierz &m1, const Macierz &m2 ){
    // ----- sprawdzenie czy m1 jest skalarem ----------------------------------
    if (m1.liczbaKolumn()==1 && m1.liczbaWierszy()==1){
        Macierz sumaMacierzy = m2;
        for(int w = 0; w<sumaMacierzy.liczbaWierszy(); w++){
            for(int k = 0; k < sumaMacierzy.liczbaKolumn(); k++){
                sumaMacierzy.at(w, k) = m2(w, k) + m1(0, 0);
            }
        }
        return sumaMacierzy;
    }
    // ----- sprawdzenie czy m2 jest skalarem ----------------------------------
    else if (m2.liczbaKolumn()==1 && m2.liczbaWierszy()==1){
        Macierz sumaMacierzy = m1;
        for(int w = 0; w<sumaMacierzy.liczbaWierszy(); w++){
            for(int k = 0; k < sumaMacierzy.liczbaKolumn(); k++){
                sumaMacierzy.at(w, k) = m1(w, k) + m2(0, 0);
            }
        }
        return sumaMacierzy;
    }
    // ----- sprawdzenie czy zgadzaja sie wymiary ------------------------------
    else if (m1.liczbaKolumn()!=m2.liczbaKolumn() || m1.liczbaWierszy() != m2.liczbaWierszy()){
        throw Wyjatek("Nie mozna dodac macierzy, niezgodne wymiary");
    }

    Macierz sumaMacierzy = m1;
    for (int w = 0; w < sumaMacierzy.liczbaWierszy(); w++) {
        for (int k = 0; k < sumaMacierzy.liczbaKolumn(); k++) {
            sumaMacierzy.at(w, k) = m1(w, k) + m2(w, k);
        }
    }
    return sumaMacierzy;
 }

//======= operator - ===========================================================
/** 
 * Funkcja przeładowanego operatora odejmowania jednoargumentowego. Funkcja zwraca
 * macierze przeciwna do macierzy podanej w argumencie - kazdy jej element jest 
 * przemnozony przez -1.
 *     @param m1  -  macierz
 *     @return -m1
 */
Macierz operator- ( const Macierz &m1){
    return m1 * (-1);
}


//======= operator - ===========================================================
/** 
 * Funkcja przeładowanego operatora odejmowania macierzy. Jezeli jedna z macierzy
 * jest skalarem to zostanie wykonane odejmowanie skalara do kazdego elementu 
 * drugiej macierzy. Jezeli zadna z macierzy nie jest skalarem i nie zgadzaja
 * sie ich wymiary funkcja wyrzuca wyjatek klasy Wyjatek.
 *     @param w1  -  odjemna 
 *     @param w2  -  odjemnik
 *     @return w1 - w2
 */
Macierz operator- ( const Macierz &m1, const Macierz &m2 ){
    // ----- sprawdzenie czy m1 jest skalarem ----------------------------------
    if (m1.liczbaKolumn()==1 && m1.liczbaWierszy()==1){
        Macierz roznicaMacierzy = m2;
        for(int w = 0; w<roznicaMacierzy.liczbaWierszy(); w++){
            for(int k = 0; k < roznicaMacierzy.liczbaKolumn(); k++){
                roznicaMacierzy.at(w, k) = m1(0, 0) - m2(w, k);
            }
        }
        return roznicaMacierzy;
    }
    // ----- sprawdzenie czy m2 jest skalarem ----------------------------------
    else if (m2.liczbaKolumn()==1 && m2.liczbaWierszy()==1){
        Macierz roznicaMacierzy = m1;
        for(int w = 0; w<roznicaMacierzy.liczbaWierszy(); w++){
            for(int k = 0; k < roznicaMacierzy.liczbaKolumn(); k++){
                roznicaMacierzy.at(w, k) = m1(w, k) - m2(0, 0);
            }
        }
        return roznicaMacierzy;
    }
    // ----- sprawdzenie czy zgadzaja sie wymiary ------------------------------
    else if (m1.liczbaKolumn()!=m2.liczbaKolumn() || m1.liczbaWierszy() != m2.liczbaWierszy()){
        throw Wyjatek("Nie mozna odjac macierzy, niezgodne wymiary");
    }
    
    // ----- wszystko ok, odejmujemy -------------------------------------------
    Macierz roznicaMacierzy = m1;
    for (int w = 0; w < roznicaMacierzy.liczbaWierszy(); w++) {
        for (int k = 0; k < roznicaMacierzy.liczbaKolumn(); k++) {
            roznicaMacierzy.at(w, k) = m1(w, k) - m2(w, k);
        }
    }
    return roznicaMacierzy;
 }

//======= operator * ===========================================================
/** 
 * Funkcja przeładowanego operatora mnozenia macierzy. Jezeli jedna z macierzy
 * jest skalarem to zostanie wykonane mnozenie skalarne. Jezeli zadna z macierzy 
 * nie jest skalarem i nie zgadzaja sie ich wymiary funkcja wyrzuca wyjatek 
 * klasy Wyjatek.
 *     @param w1  -  czynnik pierwszy 
 *     @param w2  -  czynnik drugi
 *     @return w1 * w2
 */
Macierz operator* ( const Macierz &m1, const Macierz &m2 ){
    // ----- sprawdzenie czy m1 jest skalarem ----------------------------------
    if (m1.liczbaKolumn()==1 && m1.liczbaWierszy()==1){
        Macierz iloczynMacierzy = m2;
        for(int w = 0; w<iloczynMacierzy.liczbaWierszy(); w++){
            for(int k = 0; k < iloczynMacierzy.liczbaKolumn(); k++){
                iloczynMacierzy.at(w, k) = m1(0, 0) * m2(w, k);
            }
        }
        return iloczynMacierzy;
    }
    // ----- sprawdzenie czy m2 jest skalarem ----------------------------------
    else if (m2.liczbaKolumn()==1 && m2.liczbaWierszy()==1){
        Macierz iloczynMacierzy = m1;
        for(int w = 0; w<iloczynMacierzy.liczbaWierszy(); w++){
            for(int k = 0; k < iloczynMacierzy.liczbaKolumn(); k++){
                iloczynMacierzy.at(w, k) = m1(w, k) * m2(0, 0);
            }
        }
        return iloczynMacierzy;
    }
    // ----- sprawdzenie czy zgadzaja sie wymiary ------------------------------
    else if (m1.liczbaKolumn()!=m2.liczbaWierszy()){
        throw Wyjatek("Nie mozna przemnozyc macierzy, niezgodne wymiary");
    }
    
    // ----- wszystko ok, mnozymy ----------------------------------------------
    Macierz m2Trans = m2.transponuj();
    Macierz iloczynMacierzy (vector< vector<double> >(m1.liczbaWierszy(), vector<double>(m2.liczbaKolumn(), 0.0)) );
    for (int w = 0; w < iloczynMacierzy.liczbaWierszy(); w++) {
        for (int k = 0; k < iloczynMacierzy.liczbaKolumn(); k++) {
            vector<double> wiersz = m1.wiersz(w), kolumna = m2Trans.wiersz(k);
            
            //cout<<endl<<"w  k ="<<w<<"   "<<k<<endl;
            //cout<<"wierszM1"<<w<<wiersz<<endl;
            //cout<<"wierszM2Trans"<<k<<kolumna<<endl;
            /*
            cout<<"wierszM1.begin"<<*m1.wiersz(w).begin()<<endl;
            cout<<"wierszM1.end"<<*--m1.wiersz(w).end()<<endl;
            cout<<"wierszM2Trans.begin"<<*m2Trans.wiersz(k).begin()<<endl;
            cout<<"wierszM2Trans.end"<<*--m2Trans.wiersz(k).end()<<endl;
             */
            //iloczynMacierzy.at(w, k) = inner_product(m1.wiersz(w).begin(), m1.wiersz(w).begin(), m2Trans.wiersz(k).begin(), 0);
            iloczynMacierzy.at(w, k) = inner_product(wiersz.begin(), wiersz.end(), kolumna.begin(), 0.0);
            //cout<<"iloczyn"<<w<<k<<" = "<<iloczynMacierzy(w, k)<<endl;
        }
    }
    return iloczynMacierzy;
 }

//======= operator / ===========================================================
/** 
 * Funkcja przeładowanego operatora dzielenia macierzy. Jezeli druga z macierzy
 * jest skalarem to zostanie wykonane dzielenie skalarne. Jezeli zadna z macierzy 
 * nie jest skalarem i nie zgadzaja sie ich wymiary funkcja wyrzuca wyjatek 
 * klasy Wyjatek.
 *     @param m1  -  czynnik pierwszy 
 *     @param m2  -  czynnik drugi
 *     @return m1 * inv(w2)
 */
Macierz operator/ ( const Macierz &m1, const Macierz &m2 ){
     // ----- sprawdzenie czy m2 jest skalarem ----------------------------------
    if (m2.liczbaKolumn()==1 && m2.liczbaWierszy()==1){
        Macierz ilorazMacierzy = m1;
        for(int w = 0; w<ilorazMacierzy.liczbaWierszy(); w++){
            for(int k = 0; k < ilorazMacierzy.liczbaKolumn(); k++){
                ilorazMacierzy.at(w, k) = m1(w, k) / m2(0, 0);
            }
        }
        return ilorazMacierzy;
    }
    // ----- sprawdzenie czy zgadzaja sie wymiary ------------------------------
    else if (m1.liczbaKolumn()!=m2.liczbaWierszy()){
        throw Wyjatek("Nie mozna podzielic macierzy, niezgodne wymiary");
    }
     // ----- wszystko ok, dzielimy --------------------------------------------
    Macierz ilorazMacierzy = m1 * m2.odwroc();
    return ilorazMacierzy;
}
//------------------------------------------------------------------------------
//-------------- KONIEC [funkcje pomocnicze] -----------------------------------
//------------------------------------------------------------------------------

